#include <ctime>
#include "types.h"
#include "MersenneTwister.h"

#ifndef _RANDOM_H_321_
#define _RANDOM_H_321_


static int global_seed = 5125099;
static MTRand global_rand(global_seed);

inline void reseed(int s) {
	global_rand.seed(s);
}

inline static double randdouble() {
	return global_rand.rand();
}

inline static float randfloat() {
	return (float) global_rand.rand();
}

inline static int randint() {
	return global_rand.randInt();
}

uint32 timehash(time_t t, clock_t c)
{
	// Get a uint32 from t and c
	// Better than uint32(x) in case x is floating point in [0,1]
	// Based on code by Lawrence Kirby (fred@genesis.demon.co.uk)

	static uint32 differ = 0;  // guarantee time-based seeds will change

	uint32 h1 = 0;
	unsigned char *p = (unsigned char *) &t;
	for( size_t i = 0; i < sizeof(t); ++i )
	{
		h1 *= UCHAR_MAX + 2U;
		h1 += p[i];
	}
	uint32 h2 = 0;
	p = (unsigned char *) &c;
	for( size_t j = 0; j < sizeof(c); ++j )
	{
		h2 *= UCHAR_MAX + 2U;
		h2 += p[j];
	}
	return ( h1 + differ++ ) ^ h2;
}

#define MWCMASK 0xffffffffL;

//TODO Refactor into generic abstract class that only requires randInt() to be defined

class Rand {
private:

public:
	// Re-seeding functions with same behavior as initializers
	virtual void seed( const uint32 oneSeed ) = 0;
	virtual void seed() = 0;

	// Access to 32-bit random numbers
	virtual uint32 randInt() = 0;         // integer in [0,2^32-1]
	uint32 randInt(const uint32 n);       // integer in [0,n] for n < 2^32
	double rand();                        // real number in [0,1]
	double rand( const double n );        // real number in [0,n]
	double randExc();                     // real number in [0,1)
	double randExc( const double n );     // real number in [0,n)
	double randDblExc();                  // real number in (0,1)
	double randDblExc( const double n );  // real number in (0,n)
	double operator()();                  // same as rand()

	// Access to 53-bit random numbers (capacity of IEEE double precision)
	double rand53();  // real number in [0,1)

	// Access to nonuniform random number distributions
	double randNorm( const double mean = 0.0, const double stddev = 1.0 );

};

inline uint32 Rand::randInt(const uint32 n) {
	// Find which bits are used in n
	// Optimized by Magnus Jonsson (magnus@smartelectronix.com)
	uint32 used = n;
	used |= used >> 1;
	used |= used >> 2;
	used |= used >> 4;
	used |= used >> 8;
	used |= used >> 16;

	// Draw numbers until one is found in [0,n]
	uint32 i;
	do
		i = randInt() & used;  // toss unused bits to shorten search
	while( i > n );
	return i;
}

inline double Rand::rand()
	{ return double(randInt()) * (1.0/4294967295.0); }

inline double Rand::rand( const double n )
	{ return rand() * n; }

inline double Rand::randExc()
	{ return double(randInt()) * (1.0/4294967296.0); }

inline double Rand::randExc( const double n )
	{ return randExc() * n; }

inline double Rand::randDblExc()
	{ return ( double(randInt()) + 0.5 ) * (1.0/4294967296.0); }

inline double Rand::randDblExc( const double n )
	{ return randDblExc() * n; }

inline double Rand::rand53()
{
	uint32 a = randInt() >> 5, b = randInt() >> 6;
	return ( a * 67108864.0 + b ) * (1.0/9007199254740992.0);  // by Isaku Wada
}

inline double Rand::randNorm( const double mean, const double stddev )
{
	// Return a real number from a normal (Gaussian) distribution with given
	// mean and standard deviation by polar form of Box-Muller transformation
	double x, y, r;
	do
	{
		x = 2.0 * rand() - 1.0;
		y = 2.0 * rand() - 1.0;
		r = x * x + y * y;
	}
	while ( r >= 1.0 || r == 0.0 );
	double s = sqrt( -2.0 * log(r) / r );
	return mean + x * s * stddev;
}

inline double Rand::operator()()
{
	return rand();
}

/*
 * Tausworthe Generator
 */
class TRand : public Rand {
public:
	TRand() : s(new uint32[3]), t(new uint32[3]) {
		seed();
	}

	TRand(uint32 s) : s(new uint32[3]), t(new uint32[3]) {
		seed(s);
	}

	TRand(const TRand& other) : s(new uint32[3]), t(new uint32[3]) {
		s[0] = other.s[0];
		s[1] = other.s[1];
		s[2] = other.s[2];
		t[0] = other.t[0];
		t[1] = other.t[1];
		t[2] = other.t[2];
	}

	~TRand() {
		delete[] s;
		delete[] t;
	}

	inline void seed() {
		seed((uint32)time(NULL));
	}

	inline void seed(uint32 corn) {
		float qv = sqrt(corn);
		float rv = pow(corn, 2.0);
		uint32 q = *((uint32*) &qv);
		uint32 r = *((uint32*) &rv);
		s[0] = corn;
		s[1] = ((q & 0x0000ffff) << 16) & ((q & 0xffff0000) >> 16);
		s[2] = ((r & 0x0000ffff) << 16) & ((r & 0xffff0000) >> 16);
	}

	inline uint32 randInt() {
		turn();
		return t[0]^t[1]^t[2];
	}

private:
	uint32 *s, *t;

	/*
	x_n = (s1_n ^^ s2_n ^^ s3_n)
	where,

	s1_{n+1} = (((s1_n&4294967294)<<12)^^(((s1_n<<13)^^s1_n)>>19))
	s2_{n+1} = (((s2_n&4294967288)<< 4)^^(((s2_n<< 2)^^s2_n)>>25))
	s3_{n+1} = (((s3_n&4294967280)<<17)^^(((s3_n<< 3)^^s3_n)>>11))
	*/

	void turn() {
		t[0] = (((s[0] & 0xfffffffe)<<12)^(((s[0]<<13)^s[0])>>19));
		t[1] = (((s[1] & 0xfffffff8)<< 4)^(((s[1]<< 2)^s[1])>>25));
		t[2] = (((s[2] & 0xfffffff0)<<17)^(((s[2]<< 3)^s[2])>>11));
		s[0] = t[0];
		s[1] = t[1];
		s[2] = t[2];
	}
};

/*
 * Multiply-With-Cary Generator
 */
class MWCRand : public Rand {
public:
	MWCRand() {
		seed();
	}

	MWCRand(uint32 s) {
		seed(s);
	}

	void seed() {
		x = (uint32)time(NULL) & MWCMASK;
	}

	void seed(uint32 s) {
		x = s & MWCMASK;
	}

	uint32 randInt() {
		x = (a * (x & 0xffffffffL)) + (x >> 32);
		return (uint32) x;
	}

private:
	static const uint64 a = 0xffffda61L;
	uint64 x;
};

#endif /* RANDOM.H */
